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FOREWORD 


The  work  described  in  this  report  was  performed  as  research  and  development  funded 
by  the  SLBM  Research  and  Analysis  Division  (K40)  and  the  Space  and  Surface  Systems 
Division  (K10)  of  the  Strategic  Systems  Department. 

The  development  of  this  satellite  theory  resulted  from  the  joint  efforts  of  Russell  H. 
Lyddane,  Charles  J.  Cohen,  and  the  author.  The  expertise  of  Drs.  Lyddane  and  Cohen 
strongly  affected  the  approach  taken  in  this  task,  and  their  subsequent  vigilance  provided 
a  continuous  flow  of  improvements  to  almost  all  sections  of  the  theory. 

Much  useful  information  was  extracted  from  Wayne  D.  McClain’s  two  reports 
CSC/TR-77/6010  and  CSC/TR-78/6001,  the  documentation  of  a  semianalvtic  satellite 
theory  developed  with  Paul  J.  Cefola  and  others.  We  also  had  several  discussions  with 
Dr.  Cefola  during  his  visits  to  NSWC;  his  knowledge  of  the  subject  matter  is  extensive. 
At  NSWC,  Dr.  Armido  R.  DiDonato’s  expertise  in  numerical  analysis  proved  most  helpful 
when  numerical  difficulties  arose. 

This  report  was  reviewed  by  J.  Ralph  Fallin,  Head,  Space  and  Surface  Systems 
Division. 

Approved  by: 


R.  L.  SCHMIDT,  Head 
Strategic  Systems  Department 
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INTRODUCTION 


The  K  10-Division  at  NSWC  has  many  years  of  experience  in  computing  satellite  orbits 
by  the  famed  “Cowell”  numerical  integration  methods.  For  orbits  close  to  the  Earth,  a 
stepsize  of  60  seconds  is  typical.  Increasing  accuracy  requirements  entailed  ever  larger 
gravity  fields,  and  orbit  computation  by  numerical  integration  became  more  expensive. 


In  1978,  NSWC  asked  Lyddane  to  develop  a  semianalytic  theory,  restricted  to  the  zonal 
harmonics,  but  including  Moon  and  Sun  as  perturbing  accelerations.  Lyddane  promptly 
produced  an  elegant  solution,  choosing  the  nonsingular  Poincare  elements  as  primary  vari¬ 
ables.  His  work  was  published  as  an  NSWC  TR  in  1984  (Reference  2).  As  Lyddane  points 
out,  such  a  solution  should  be  at  least  two  orders  of  magnitude  faster  than  the  customary 
Cowell  integration. 

Subsequently  the  decision  was  made  to  develop  a  more  complete  semianalytic  theory 
here  at  NSWC.  Lyddane  and  Cohen  offered  their  support,  a  most  welcome  event.  The 
major  feature  of  the  mathematical  model  was  the  addition  of  the  tesseral  gravity  field, 
with  particular  attention  given  to  the  notoriously  troublesome  resonance  perturbations. 
This  phase  of  the  work  has  now  been  completed,  and  the  computer  program  is  undergoing 
checkout. 


In  assembling  the  program,  the  tesseral  algorithm  was  merged  with  Lyddane's  op¬ 
erational  zonal  routine.  Since  the  latter  is  documented  in  Reference  2,  this  report  will 
describe  the  effects  of  the  tesseral  field  and  explain  how  various  types  of  terms,  especially 
those  with  small  divisors,  are  treated. 


TESSERAL  HARMONICS  IN  POINCARE  VARIABLES 

THE  DISTURBING  FUNCTION 

In  this  section,  the  disturbing  function  for  the  tesseral  gravity  field  will  be  derived  in 
terms  of  slightly  modified  Poincare  orbital  elements. 
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The  Earth’s  potential,  as  adopted  by  the  International  Astronomical  Union  (Refer 
ence  1),  is  written 


t'  =  ^  +  £E(v)"p»  m  (sin  f3)  [Cnm  cos  m\  +  Snm  sin  mA]  1 

r  ^  „=2  m= o  '  7  '  ' 


For  the  sake  of  brevity  in  subsequent  derivations,  we  write  the  potential  in  the  form 


Figure  1.  Oibital  geometry 
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equations  and  figure  above, 

P 

=  k2M 

k 

—  Gaussian  constant 

m 

=  mass  of  Earth 

Re 

=  equatorial  radius  of  Earth 

n 

=  degree  of  tes serai  harmonic 

m 

=  order  of  tesseral  harmonic 

Pntn 

=  associated  Legendre  function 

A 

=  longitude  from  Greenwich 

0 

=  latitude  (declination) 

9 

=  colatitude 

Cnmi  Snm 

,=  unnormalized  gravity  coefficients 

7nmi 7 nm 

=  complex  form  of  C,  S  and  its  conjugate 

i 

=  n/-1 

9 

=  inertial  rotation  rate  of  Earth 

et 

=  GMST  =  Greenwich  Mean  Sidereal  Time 

A  comparison  of  (1)  and  (2)  yields  the  relations 


Cnm  —  7nm  "4“  Tr 


Snm  —  *(7nm  7  nm) 


CnO  —  7«0  T  TnO  — ' 


7nm  —  2(^nm  *  Snm)  1 


7  nm  —  o  (^nm  d~  *  Snm) 


(3) 


Let  the  Hamiltonian  be 


F  =  F0  +  R 


where 


F0  = 


(4) 


( 5n  | 


*=-EE(t)  P„m(cos«)  [i„„,e""A  +  =f„, 

r  n=2  m  =  0  '  ' 


—  imAl 


(oh) 
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Since  the  disturbing  function  R  contains  the  time,  the  Hamiltonian  F  is  no  longer 
a  constant.  This  feature  does  not  affect  the  validity  of  our  theory.  The  Delaunay  and 
Poincare  variable  L  is  defined  by 


L  -  y/jia 


where 


a  =  semimajor  axis 


In  preparation  for  subsequent  developments,  the  disturbing  function  will  be  clad  in  a 
differenct  form.  First  we  note  that 


Pnm(x)  =  (1  —  a*2)m/2 


dmPn(x ) 
dxm 


and  introduce 


Pjjm{x) 


dmPn(x) 

dxm 


In  our  application,  this  yields 


Pnm(cos6)  =  (sin  6)m Pnm(cos  6) 


(6) 


Define 


T„m(0,  A)  =  P„m(cos0)e,n,A 

=  Pnm(cos0)  [sin0e'Ajm 
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Also  substitute 


a 

r 


,n+ 2 


L2n+2 


Rr 


(8) 


Finally,  let 


Rnm  — 


,n+ 2 


L2n+2 


/  n  \  n+i 

RnE{-)  Tnm(8,  A) 


(9) 


With  these  modifications,  the  disturbing  function  becomes 


n 

R  ^  ^  ^  ^  [inm  Unm 

n=2  m=0 


d  T  nrn 


Since  the  zonal  harmonics  have  already  been  done  by  Lyddane  (Reference  2).  we  will 
separate  R  accordingly  and  write 


R  '  ^  '  CnO-^uO  d  ^  ^  ^  [7nm^mn  d  Tnm^i1”] 


n=2 


n=2  m  =  l 


(10) 


Next  we  wish  to  express  the  surface  harmonics  in  terms  of  the  orbital  elements 
u  —  f  +  g,  h,  and  I .  They  are 


u  =  argument  of  latitude 
/  =  true  anomaly 
h  —  longitude  of  ascending  node 
I  =  inclination 
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In  doing  so,  we  follow  the  work  of  Allan  (Reference  3)  who  appeals  to  Izsak  (Refer¬ 
ence  4).  Thus 


Tnm(9,  A)  =  J2  Fn mp(I)  e«[( n-2p)«+m(ft-*0] 

p=0 

with  the  inclination  function 


Fnmp(I)  =  I"' 


(n  -f-  m)! 

2  np!(n  —  p)! 


n=k  i 


c3n  —  m  —  2p—2kgTn  —  n  +  2p+2k 


where 


=  A(n  —  m  —  2p-fjn  —  m  —  2p|) 
A;2  =  |(3n  —  m  —  2p  —  \n  +  m  —  2p|) 


I 

c  =  cos  — 


(11) 


Let  us  introduce 


v  =  n  —  m  —  2p 


and  a  slightly  modified  inclination  function 


(n  +  m)!  /  n  aY2n  ~  2/A  ( 

2 np!(n  —  p)!  \  k  J  \n 

k=k\ 


2  P  \  c2n  +  v-2k  ,,-v+ik 

rn  —  k  J 


(12) 
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where 


h  =  \(v  +  M) 

k2  -  |(2 n  +  v  -  |2 m  +  v\) 
Then  Tnm  may  be  written 


T.m(«,  A)  =  in'm  ^  F-mp(c, s)  e<[(»-2p)-+-»(‘-«)l 

P= 0 

or,  putting 

jn  —  m  _  ei(n-m)ir/2 


Tnvi(e,\)  =  YJFnmp(c,s)e'^  (13) 

p=0 

where 


V’nmp  =  (n  -  2p)u  +  mh  -  mOt  +  (n  -  in)-  (14) 

Next,  working  toward  Poincare  variables,  we  express  t/>nmp  in  terms  of  h  and  u  4-  h 
for  direct  orbits  and  of  h  and  u  —  h  for  retrograde  orbits.  Accordingly,  we  write 


Tnm  =  Tnmp 

p=0 
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where 


Tnmp 


—  Fnmp^i  ■s) 


ei[(n—  2p)(u+h)  —  (n  —  m—  2p)h  —  mOt  +  (n  —  m)n/2] 


(15) 


for  direct  orbits  and 


T  —  r  (n  J[(n-2p)(u-h)+(n+m~2p)h-mdt+(n-m)n/2] 

•* nmp  —  *  nmp\ci  c 


(16) 


for  retrograde  orbits. 

For  direct  equatorial  orbits,  s  =  0  and  c  —  1 ;  for  retrograde  equatorial  orbits,  the 
values  interchange,  becoming  s  =  1  and  c  —  0.  One  would  therefore  look  for  symmetry  in 
the  exponents  of  c  and  s  to  overcome  the  indeterminancy  of  h  in  equatoxial  orbits,  both 
when  c  =  0  and  s  =  0.  It  turns  out  the  symmetry  is  there  if  Fnmp  is  compared  with 
Fnm,n-p  and  then  Rnmp  with  Thus  for  Fnm<n-p,  if  we  replace  the  index  of 

summation  k  by  k'  =  n  —  m  —  k ,  we  may  write 


2  (n  -p)\p\  kf^k,  \n 


2  P 

in  —  k' 


2n  2p\  -»+2k'  ,2n+v-2k' 

k'  1 


=  (-l)n  m Fnmp{s,  c) 


(17) 


since 


k\  =  n  —  m  —  k2{n,m,n  —  p)  =  k\(n,m,p) 


and 


k'2  =  n  —  m  —  k{  ( ??.,  m,n  —  //)  =  A:2 ( 77 •  m,p) 
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For  the  retrograde  orbit  now  consider  Tnm,n-p  •  From  (16),  we  have 


rp  _p  (  \  i[-(n-2p)(u-h)-(n-2:p)()l-vh  +  (.n-m-2p)et+(n-m)w/2\ 

*  nm,n  —  p  —  1  nm,n— pvc»  d )  c 

_  (—l')n~mFnmp(s  c)  e-l!(n_2p)(u_/l)+l/,l+m®t]gi(n-m)7r/2 


after  using  (17).  Since 


(_l)n-,n 


ei(n-m)ir/2 


—  e-'(n_m)7r/2 


we  obtain 


T 

nmyn—p 


nmp 


(s,c)e 


—  i[(n  —  2p)(u—h)  +  vh-\-m8t+(.  n  —  m)n/2] 


This  form  for  retrograde  orbits  is  conveniently  close  to  the  form  for  direct  orbits,  (15), 
and  henceforth  (18),  together  with 


T 

-*■  nr 


=  ET' 

p=0 


nm,n-pi 


(19) 


will  be  used  for  retrograde  orbits.  We  will  return  to  the  effects  of  summing  on  p  vs  n  —  p 
when  we  later  compare  Rnmpq  and  Rnm  n-Pt-q ,  in  (39). 

Now  let 


e 


+  1  for  I  <  y± 
-1  for  /  >  f± 


(direct  orbits) 
(retrograde  orbits) 
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and  also 


h'  -  eh 


H'  =  eH 


c 


i 


(21) 


Let’s  also  define  a  new  function 


V’nmp 


e 


—  2 p'j  (u  +  h'  —  eOt'j  —  (n 


m  — 


or 


^nmp  e 


fn  -  2 pj  (n  -f  h'^j  -  vti  -  emQt  +  (n  -  rnj  ^ 


(22) 


With  these  preparations,  we  obtain  a  single  equation  for  Tnm  good  in  both  the  direct 
and  retrograde  cases: 


Tnr 


-E* 

P= 0 


nmp 


(c',s')e 


tlb> 

t  nr 


(23) 


Let’s  verify  this  assertion: 


€_—  ■ +1 


-Eft 

p=0 


mp 


(c  s)  e'[{n-2p){u+h)-''h-mkMn~m)^ 


which  does  agree  with  15. 


10 


NSWC  TR  89-109 


£  =  -1 


^  ^  Fnmp  ($i  c)  e 
P= 0 


t[(n— 2p)(u  —  h)  +  vh  +  m9t+(,n  —  m)i 


which  is,  indeed,  (18). 

Now  recall  i?n„, ,  given  by  (9): 


Rnm  — 


n+2 


£2n+2 


/  a  \  1 

RnE[~)  Tnm{6,\) 


This  becomes 


R 


nm 


~^RnE  {^Y+1  ^2FnmP(c\s')  ei^-2p)(g+h')-,h' -tmet+(n-,n)  f]  ei<(n-2p)f 
P~  0 


From  this  expression  we  extract  the  short-periodic  content  and  expand  the  latter  in  a 
Fourier  series  in  the  mean  anomaly  l  with  Hansen  coefficients  X  (/?): 


i  ,  w 

^ )"  eif{n~2p)f  =  Xn-2pfq~2v  (P)  eu{"-2p+q)l 

—  OO 


(25) 


where 


0  = 


e 

1  +  \/l  ~  e2 


(26) 


The  e  in  (26),  as  well  as  below,  is  the  eccentricity.  Since  the  base  of  the  natural  log 
always  involves  i ,  confusion  should  not  arise. 

Also  note  that  the  validity  of  (25)  for  e  =  —  1  is  easily  established  by  taking  it  to  be 
the  complex  conjugate  of  the  f  =  +1  case. 
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If  we  write 


Rnm  —  'y  ]  y  ^ 


p= 0  q-=  —  oo 


we  obtain 


_  if _ pn  p  U[(n-2p)(g+h')-vh'-(m8t+(n-m)Z  ] 

~  £2n+2  IXErnmp  Vc  t  *  J  e 

*  X„-:*r2P  («  e1'<"-2',+*M 

n-f  2 


where 


=  €  (n  -  2p  +  g)  /  +  (n  -  2p)  (fir  +  /»')  -  vti  -  emdt  -f  (n  -  m)  £  (28a) 


In  preparation  for  introducing  the  Poincare  variables,  rearranging  terms  yields 


KmPq  =  e  (n-2p  +  q)(l  +  g  +  h')-q{g  +  ti)-vti  -emW  +  (n-m)J  (29) 


Our  Poincare  variables,  modified  to  accommodate  all  inclinations,  are  defined  by 


L  A '  =  l  +  g  +  h 

£  =  \/2.{L  —  G)cos(<?  +  h')  q  =  —  v/2 (L  —  G)  sin (g  +  h') 


a  =  v/2(G  -  H')costi  t  =  v/2(G  -  H')smh' 
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Also  anticipating  further  developments,  we  define  the  functions 


Fnmp  (c  5  s  )  —  rjr Fnmp  (c  )  s  ) 
(s')1 


Then,  in  (28), 


F„mr  (c\  s')x:^X'2p  W  tiK~’'  =  (s')w  A™,  (c',  s’)  I3i’'x;:;n'£-2p  w 

+  eu[(n-2p+?)A' -q(g+h')-vlt' -(in«(  +  (n-m)f] 

=  F„mr  (c\s')  W  (*')m 

*  e~ituh ) 

+  e»<[(n— 2p+q)X'  —  tm6t+(n—  m)y] 

*  fie~i({9+h')!,9n<l  ^ 

ie[(n-2p+g)A'-fmi9<  +  (n  —  »«)■§•] 


For  the  remaining  steps,  as  well  as  for  the  computer  program  algorithm,  the  following 
auxiliary  relations  will  be  useful; 


2{L  —  G)  =  (2  +  j/2 


2(G  -  H')  =  cr2  +  t2 


G  =  L  -  i  (e  +  .)2) 
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0  = 


s  = 


( 2  +  *i 2 


4 I-(£2+r?2) 


<72  +  t2 


4L-2(e+V2)) 


Then 

s'e^'^^s'ie-^Y39^ 


=  5  [cos  h  —  t  sm  h  j 
=  s'  [cos h'  —  ie sgn  v  sin h'] 


/  *2+r2 

o  +  i  esgnv  t 

V  4L-2(£2  +g2) 

sjo2  +  r2 

juk'sgn *  =  <r  +  (tesgnu)T 

y/4L-2(?  +  if) 


Quite  similarly  we  obtain 


_  Z  +  (iesgnq)n 
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Substituting  (36)  and  (37),  we  reach  teh  following  form  for  i?nmp9 ,  free  of  small 
divisors: 


,n+2 


R 


nmpq  £2n+2 


R%R nmp  (c\  S  ) 


<7  +  i  er  sgn  v 

[v/4L-2(e+.F)J 


£  +  iei]  sgn  q 

W+W) 


i  l?l 


3t« [(n— 2p+q) A'  +  (n  —  m)  ^  ^  — 


(38) 


We  note  that  Rnm  (c  =  —1)  =  Rnm  (e  =  +1)  but,  because  of  a  shift  in  the  summation 
indices  p  and  q,  we  have  Rnmpq  (e  =  —1)  ^  Rnmpq  (e  =  +1)-  To  see  the  effect  of  the  shift, 
we  return  to  (28): 


Rnm „  (e - 1)  =  R%Fnmp  (s,  c)  W 

*  e-i[(n~'iP+,l)l+(n-‘iP')(-9-h)+i'h  +  rn8t+(n-m)l] 


n+2 


Rnm,n-p,-q  (e  —  +1)  —  £2n+2  ^E^nm,n-p  (s,  c)  -‘^_(n-2p+j)  P  (/^) 

^  et[(  —  u+2p— $)/+(  — n+2p)($+ft)  +  (n  +  m-2p)/t  —  m0f  +  (n-m)^] 


n+2 


V — n  — 1 ,  — (n  — 2p) 


using  (17) 


/d‘+2 

£2n  +  2 


— (n— 2p+fl)/— (n— 2p)(j  — ft )  — I'ft  —  m^<  +  (  n  — m )  ■£] 
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Using  (-I)"-”*  =  ei(n-m)n 


n+2 


*  et[(n-2p+9)/+(n-2p)(fl-fc)+x/fe+met-(n-m)fjei2(n-m)f 


implies 


R 


nmpq 


(e  =  -l)  =  R 


nm,n— p, 


(e  —  +1) 


(39) 


Note  that  use  was  made  of 


y  —  n—  1  ,n— 2p  y  — n  — 1,  — (n— 2p) 

^n-2p+,  -^_(„_2p+g) 


which  follows  from  (25)  and  its  complex  conjugate. 


SUMMARY  OF  RESULTS 


We  recapitulate  showing  the  Hamiltonian  framework  and  the  definitions  of  the  vari¬ 
ables  that  enter  Rnmpq  •  The  canonical  variables  are: 


L  A '  =  l  +  g  +  h 

£  =  \Z^{L  —  G ) cos ( g  -f  h')  t)  =  -y/2 (L  -  G')sin(g  +  h ')  > 

a  =  y/2 (G  -  H')costi  r  =  -y/2(G  -  H')smh' 


(40) 
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Note  that 


H'  =eH  ti  =  eh 


where 


e  = 


+1  for/<f± 


-1  for  I  >  f  ± 


Once  set,  e  is  never  changed  during  a  run.  The  Hamiltonian  is 


F  =  F0  +  R 


where 


F0  = 


jF_ 

2  L2 


Mil  Rnm  4~  rnRnm] 


n= 2 


n=2  m=l 


n  oo 


Rnm  ^  ^  ^  ^  Rnmpq 

;>=0  q=  —  oo 


„"  +  2  _ 

n  _  _ pn  T?J  I  J  J\ 

nnmpq  -  L2n+2  Er  nmp  ’  *  / 


a  +  (i  e  sgn  v)  t 

v/4i-2u2  +  //:fy 


(41) 


(42a) 


(42b) 


y-n-l,n-2/>  ,  ,, 
n  — 2;>+ij  ON 


£  +  ( iesgnq)i / 

\/4L-(e+>)2) 


l?l 


,i<[(n-2p+<f)A'+(n-»ii)f  ]  -i»n0« 


(43) 
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Note  that 


v/4L-2(e+*?2)  =  2v/G 


y/4  L-(?+r,*)  =  y/2(L  +  G) 


FJ 

A  nmp 


(c\*')  =  (~ 1)6 


(n  -  r')i(n  +  r')!  pa,6  , 

2n(n  —  m)!p!(n  —  />)'.  j  "-r'1  j 


where 


r  =  n  —  2p 


a  =  m  —  r 


a  =  |m  —  r| 


6  =  ^(a'  -  a) 

r'  =  ^(a  +  &) 

/  /,  g2  +  r? 

V  4Z  —  2(£2  4-  rj>2) 


fe  =  |m  +  r|  A"'  —  2(c')2  —  1 


Compute 


p0aV)  =  i 

p^(x')  =  -(/?  +  l)x'  +  -(a  -  b ) 


where 


p  —  (i  ~h  b  -f~  1 


(44) 
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V  (*•)  -  2d{d  +  p  -  l)(2d  +  p  -  3)  {  ^  +  "  -  2) 

*  [(2d  +  p  -  3)  (2d  +  p  -  1)  x'  +  (a2  -  62)] Pad'\  (s') 

-  2(d  +  a  -  1)  (d  +  b  -  1)  (2d  +  p  -  1)  Pad'\  (*')}  (45) 


Starting  with  d  —  2 ,  continue  until  d  =  n  —  r' . 


y  —  n  —  1  ,n  —  2p 
— 2p+g 


(/))  =  (-i)“ 


(i  +  /?T 

(1  -  /?2)2n-l 


Y,MsNt(12j 
}= 0 


(46) 


where 


/?  = 


£2  4.  ??2 


4L-K2  +  ,2) 


e 

1  +  >/l  -  e2 


=(-ir'r(“ 


M3  =  (-1) 


—  tA 


2-— '  —  u  /  v! 

r=0 


5  =  j  +  ^(M  +  9) 
t=j  +  l(\q\-q) 

a  =  2n  —  2p  +  q  —  1 
7  =  2p  —  q  —  1 
b  =  (n  -  2/;  4-  <?)  \/l" -  e2  , 


(47) 


(48) 
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Alternatively,  Ms  and  Nt  may  be  computed  using  Lyddane’s  relations: 
Starting  with 

M  (0,  a;  b)  =  1  'j 

M  (s,s  —  1;6)  =  (— 1)3-t  | 

s!  ' 

compute 

Ms  =  M  (s,  a;  b)  =  M  (s,  a  —  1;  b)  —  M  (s  —  1,  a  —  1;  b) 

Nt  =  N(t,r,  b)  =  M(t,r,-b) 

The  last  line  means  that  N  is  obtained  similarly  to  M,  but  writing  -b  for  b. 
To  truncate  the  summation  in  (46),  let 

Aj  =  \M3Nt\p> 

Terminate  when 

Aj  +  Aj+ 1  <  max  ( A* )  •  tx ,  k  <  j , 


where  tx  is  an  input. 
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EQUATIONS  OF  MOTION 


THE  DERIVATION 

The  basic  equations  of  motion  for  our  canorical  Poincare  variables  take  the  well-known 

form 


OF 

{  1 

dF 

d\' 

A  = 

dL 

dF 

dF 

dr) 

7)  = 

~  dt 

dF 

dF 

dr 

T  = 

da 

(54) 


where  F  is  the  Hmailtonian,  developed  in  the  previous  section.  In  this  section  the  partial 
derivatives  indicated  in  (54)  will  be  derived. 

As  mentioned  before,  the  central  force  field  and  the  zonal  harmonics 


h^c 

71=2 


n0 


R»io 


have  already  been  treated  by  Lyddane  (Reference  2).  Hence,  we  will  address  here  the 
tesseral  part 


AF27’  =  R2T1 


(55) 


where  the  subscript  2  designates  that  all  terms  are  of  order  two  or  higher.  Following  (42b). 
we  write 


R'2T  —  ^  ^  ^  [7um-Riu 


T  7 11  ui-^  mu 


n= 2  m=l 


(5C) 
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Also  from  the  previous  section,  namely  (43),  each  individual  term  in  (56)  is 


when 


,,»+2  _ 

ft  __  _  pJ  (J  J\ 

■^nrnpq  jJZn+2  r  nnip  Vc  i  *  ) 


a  +  (i  e  sgn  v)  t 

v4L  — 2  0e2+>F) 


n_2p  J. 

yin-2p+q 


j  +  (iesgnq)i) 
y/*L  -  W  +  *?2) 


l?l 


0t't[(n— 2p+g)V+(n— m)-|]  -im8t 


(57) 


(•s')2  = 


/?  = 


<t2  +  r2 


4L-2(£2  +  >/2) 


=  ^l-ls')2 


£2  +  if 


AL~2(e  +V2) 


(58) 


Introducing  short-hand  notation  for  (57),  we  write 


R  =  .\FB'^XCWD 


(59) 


where 


A  =  A(L)  X  —  X(L,  £,  ??) ' 

F  =  F(L,^,7/,a,r)  C  =  C(L,£,»>)  * 


B  =  B{L,  <7,  r) 


D  =  D(X'J)  ) 


(60) 


Now  let  y  designate  any  one  of  our  Poincare  variables.  Then  eventually 


OR  _  RdA  RdF  RdB  RdX  R  dC  ROD 
dy  ~  A  dy  +  F  dy  +  dy  +  x  ~dy  +  l<?lC  +  D  dy 
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Hence,  we  need  the  partials  of  the  factors  in  (60)  with  respect  to  their 
arguments.  In  above  sequence, 


L 2"+2 


which  yields 


dA 

dL 


=  -(2n  +  2)j 


Next 


F  =-(-!)' 


(h  —  ?'')!(n  +  r')! 
2n(n  —  m)!p!(n  —  p)! 


(cT  K-r'  (*') 


r  = 

—  2p 

6 

=  \(a'  -  a) 

a'  ~ 

m  —  r 

r' 

=  + 

a  = 

|  m  —  r| 

P 

=  a  -t-  b  -j-  1 

b  = 

|m  +  r| 

x1 

=  2(c')2  —  1 

)2  = 

(T2 

41  -2 

+  r2 

(t2  +  V2) 

c' 

Also  note 


4L  -  2(£2  +p2)  =  4G' 


Poincare 


(Gla) 


(62) 


(63) 


(64) 


(64a) 
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We  also  need 


n*'V)  =  i  (ss) 

Pf'V)  =  i(?+l)*'  +  i  (<.-!>) 

<*'>  =  { (M + '  ~ 2) 

*  [(2d  +  p-3)(2d  +  p-  l)x'  +  (a2  -62)]P;i61(x') 

-2(d  +  o-  l)(d  +  6-  l)(2d  +  p-  1)P£.*2  (*')}  (66) 


For  convenience,  let’s  temporarily  drop  a  and  b  from  P,  the  prime  from  c',  s' ,  and 
x1 ,  and  define 


A'  =  (  11*  (tt-r')!(n  +  r')! 

2"  (n  —  m)\p\  (n  —  p)\ 


Then 


F  =  KcbPn-r'  (x) 


(67) 


implies 


dF  _  d£dc 

dc  dy 


(68) 
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where  y  is  any  one  of  our  six  Poincare  elements.  Using  (67),  first  do 


~  =  Kbcb-lPn_r,  (x)  +  I<cbdP-n-r 
dc  ax  dc 


(69) 


But  since 


x  ~  2 c2 


-1, 


and  (69)  becomes 


_  (  x)s  ( n  -  r')i(n +  r')i 
dc  2  "(n  —  m)!p!(n  —  p)! 


6Pn_r.(x')  +  4c2 


dP n—  r'  (*£  ) 

dx 


(70) 


In  order  to  obtain  ,  write  (66)  in  the  abbreviated  form 


Pd(x)  =  —  {a2[a3x  4-  a4]Pd_1(x)  -  a5Pd_2(x)} 


implies 


dPd{x) 

dx 


—  {  a2  [a3x  +  a4] 

«i  l 


dPd~\{x)  dPd_2(x)' 1  a2a3 


dx 


-  a  5 


dx 


+  —Pd-M 

ai 
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or 


dPd(x) 

dx 


_ 1 _ 

2d  (d  +  p  —  1)  (2d  +  p  —  3) 


{ 


(2d  +  p-  2) 


*  [  (2d  +  p  -  3)  (2d  +  p  -  1)  x  +  (a2  -  62)] 


-2(d  +  a-l)(d+b-l)(2d  +  p-l)  dPd'2^ 

dx 

(2d  +  p  -  l)(2d  +  p  -  2)  \ 

+  — 2d(d  +  P-i) — Pi~'{x)S 


(71) 


From  (65)  we  derive 


dP0(x) 

dx 


dPx(x) 

dx 


5('>  +  1) 


(72) 


Note  that  (71)  has  the  form  of  (66),  with  P  replaced  by  ,  plus  an  additional  term. 
The  recurrence  loops  for  P  and  ^  should  probably  be  computed  simultaneously. 

Now  turn  to 


dc 

dy‘ 


Since  we  have  now  introduced  s  in  addition  to  c,  in  (64),  let’s  make  use  of 


dc 

dy 


-1  d_ 

2c  dy 


(73) 
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The  five  partial  drivatives  are  found  to  be 


dc 

s2 

dc 

-fc2 

dL 

~  2 cG 

II 

1 

l<o 

2  cG 

dc 

—a 

dc 

—  1)S 2 

da 

~  4  cG 

dt) 

2  cG 

dc 

—  T 

dr 

~  4cG 

(74) 


The  final  results  for  F  are: 


dF 

~dl 

dF 


{s')2  dF 
2d  G  dd 

(s')2  dF 


d£  ^  ®  2 dG  dd 

°L  =  (-n)(j!t£ 

dr)  y  ”  2 dG  dd 


dF  _  -a  dF 
da  ~  4dGdd 


d£  _  -r  dF 
dr  4c' G  dd 


(75) 
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where 


OF 

dd 


(-1) 


s  (n  -  r')\  (n  +  r')\ 

2n  (n  —  m)lp\  (n  —  p)\ 


V) 


6-1 


* 


dP‘t,  (X1) 


ip;i‘r.(x')  +  4(c') 


dx' 


(76) 


with 


Pn-r‘  (x>)  defined  by  (65)  and  (66), 


dK  iV  (*’) 


dx' 


defined  by  (71)  and  (72)  . 


The  next  factor  is  B .  With  later  developments  in  mind,  write  it  as  follows: 


a  +  (i  e  sgn  u)  r 

y/4L-2U*  +  n1) 

cr  +  (ie  sgn  v)  r 

27S 


28 


NSWC  TR  89-109 


From  (77),  we  immediately  obtain 


dB  _  1 

~da  ~  2s/G 

dB  ( iesgnu ) 

d7  “  2  VG 


(78) 


For  the  other  three,  write 


B  =  -  [a  +  (ie  sgn  v)  r] 


g-i/2 


Hence 


dB_  _  5B5G 

5G  5y 


First  get 


or 


dB  _  1  [a  +  ( i  e  sgn  u)  r] 

dG  ~  ~  2  2\[GG 

dB  _  -B 
dG  ~  26' 


(SO) 


and  (79)  becomes 


dB  _  -B  dG 
dy  ~  2 G  dy 


(SI) 
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Since 


G=L-\(e+j) 


dG 

dL 


=  1 


dG 

dt 


=  -e 


dG 

dr) 


= 


Therefore 


dB_  _  -B 
dL  ~  2 G 

dB_ 

di  ~  2  G 

dB__rjB_ 
drj  ~  2G  . 


(82) 


Let’s  turn  to 


X  =  X(L,t,v) 


Using  (46)  through  (48),  we  write 


X  =  E  MsNt/32j 
j= o 


(S3) 
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where 


E  =  (-i)i,i  (1±£T_ 

(1  - 


1  +  v/l3"? 


P2  = 


e  +  n2 

4  !-(£*  + r/2) 


M,  =  (-iy±(a-v)b^ 

'  \S  —  V  J  v\ 

u=0  x  7 


*  =  ;  +  + 


*  =  j  +  9  (M  ~9) 


a  =  2n  —  2p  +  q  —  1 
7  =  2p  -  9  -  1 


6  =  (n  -  2p  +  q)  1  -  ^  (£2  +  t?2  ) 


=  (n  -  2/j  +  g)  \/l  -  e2 
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Writing  termporarily 


Q  =  n  -  2p  +  q, 


other  useful  relations  are 


b  = 


Q 


1-p2 

l  +  P2 


2/? 

6  l  +  P2 


Now  we  proceed  to  obtain 


dX  dX  dp2  r  , 

dy  ~  dp2  dy'  V~  '  ^  V 


(87) 


(88) 


Starting  with  (83): 


dX 

dST 


uE  ^  w  (dMs  „  w  dNt 

~  dP2  ^  Ms  t/3  +E^0  \~dP*Nt  +Ms'dP2 


+  EYJM,Nt]{p2) 

i= o 


2\J_1 


(p2Y 


(89) 


The  derivatives  in  (89)  are  as  follows.  First 


dE 

dp2 


n 


+ 


2  n  -  1 


(l  +  P2)  (1  -P2) 


E 


(90) 
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Next 


dMs  _  v\  v^V  1  db 

dp2  —  v )  i?!  dp2 

dNt  ,  ( 7  -  w\  w(-b)w-\- 1)  db 

d/32  1  j  ^  W-uJ  u>!  dp2 

w=  1  N  ' 


(91) 


Since 


b  =  Q 


\-p2 

1  +  P2 


db  _  -2 Q 

dp2  ~  (1  -p2)2 


(92) 


(91)  becomes 


dMa 

dp2 


-2  Q 


(i  +  P2Y 


(-d'E 


r=l 


a  —  v 
s  —  V 


b^v 

b 


dNt  _  +2Q  t  fr-wy-br  rv 

dP2  (1  +  P2)2  "  \f-uv  u>!  (-6) 


or,  since 


Q  =  l  +  /?2 

6  1  -  p2  ’ 
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dM3 

Up 


dNt 

dp 


-2 

l^P 


(-d*E 

U  =  1 


1  -  /)'  ’  ^{t-wj 

W=l  x 


■w 


(93) 


Substitution  of  (90)  and  (93)  into  (89)  yields 


dX 

dp 


n  2n  —  1 

(1  +  P)  +  l^P\ 


(_1)M  Ji±^!L  TMaNtpi 
(1  ~P)  u 


+  (-l) 


!?l 


(1  +  P) 


2  \  n  °° 


(1  -P) 


2n-l  Z_^ 


7=0 


-2 


[l-p 


(-1) 


s  —  V  )  v\ 


tu=l  x 


U) 


P1 


+(-x)h' a1-;?- 


(04) 


Let  now 


(M„)  =  (-!)*£ 

U=1 


t 

{Nw)  =  (-\)tYJ 

w=  1 


7  -  «A  (-6) 


t  —  w 


w 


W ! 


(95) 


34 


NSWC  TR  89-109 


Then 


dX 
dd 2 


+ 


2  n  -  1 


(1+/92)  '  1-/92 


n  —  1 


X  -2(-l)l?l 


(1+^2) 

a-^2)2" 


Nw) 


+  (Mv)Nt]02j  +(-l)l?l 


(1  +  /52)" 

(l_^2)2n-l 


hXtJ32}~2 

j= i 


(96) 


Next  we  differentiate 


4L-e-  rj2 

with  respect  to  £,  and  It  will  be  useful  to  note  that 


(52  +  >)2)  =  2  (I-G) 
41 -({2  +  t;2)  =  2(1  +  0 


Then 


di92  2{G  -  L) 
dL  ~  {L  +  G) 2 

5/92  2  £  I 

aT  ~~  (i  +  G)2 

chi2  2i]  L 

d>i  ~  (L  +  G)2  . 


(97) 


35 


NSWC  TR  89-109 


Using  (SS)  and  (97),  we  finally  obtain 


OX  2{G-L)dX\ 


dL  ' 

~  (L  +  G)2 

a/?2 

dX 

2  iL 

ax 

dZ  = 

'  ( L  +  G )2 

a/?2 

dX  _ 

2  7/1 

ax 

di] 

"  (T  +  G)2 

a/?2 

(98) 


with  dX  /  d02  given  by  (96). 

An  alternative  path  to  obtain  the  partials  given  in  (98)  uses  the  recurrence  relations 
proposed  by  Lyddane.  Instead  of  (84),  which  uses  potentially  troublesome  binomial  coef¬ 
ficients,  we  write 


Ma  =  M  {s,  a;  b)  =  M  (s,  a  -  1;  b)  -  M  {s  -  1,  a  -  1;  6)' 
Nt  =  N  (t,  7;  b)  =  M  (t,  7;  -6) 


(99) 


The  recurrence  (99)  uses  the  starting  values 


M  (0,  q;  6)  =  1 


M  (s,  s  -  1;  6)  =  (-1)3  — 

s! 


(100) 


Equations  (99)  state  that  Nt  may  be  obtained  using  the  algorithm  for  A/,  by  substituting 
t  and  7  for  s  and  c»  and  by  writing  -b  in  place  of  +6. 
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Dropping  a  few  subscripts,  we  copy  from  Lyddane 


dM  1  i  ,x 

—  =  -M  (a  -  1,  a  -  1;  b ) 

dN 

—  =  +M(<-  1,7-1;  -&) 


implies 


dM  dM  db  +2Q 


d/32 

db  dp2  (1  +  p2j 

diV 

-2Q 

dp2 

"  (l  +  /?2): 

M  (s  —  1,  a  —  1;  6) 


M(t  -  1,  7  -  1;  -6) 


We  now  return  to  (89)  and  substitute  (90)  and  (101): 


dX 


n 


+ 


2  n  -  1 


dp  i(i+ n  (i-P) 


X  +  2(-l)M(n-2p  +  q) 


(l+P) 


n—2 


(1 


*  M  ( t ,  7;  —b)  —  M  (s,  a;  b)  M  (t  —  1,  7  —  1;  -6)]  f32j 

(l  4.  /92\n  0= 

+  (-1)'"  «■:  7;  -H 

This  equation  takes  the  place  of  (96)  in  (98). 

The  fifth  factor  in  the  disturbing  function  is 


-f  ( iesgnq)i / 

v/4I-(e  +  r;2) 


(101) 


l,a-l;  6) 


(102) 


(103) 


37 


SUMMARY  OF  RESULTS 

Throughout  this  section,  the  abbreviation  R  was  used  for  the  individual  term  Rn  pq 
As  in  (59),  we  write 


R  =  AFB^X  C^D  (106) 
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where 


]in+2RnE 

~  £2n+2 


=  -4  (i-i 


F  =  (-!)' 


(n  —  ?*')!  (n  -j-  r/)? 

2"  (n  —  m)!p!  (n  —  p)! 


{x')  =  F(L,  £,  rj,  a,  t) 


<r  +  (i  e  sgn  v)  r 


X 


(_1}M 


(l  +  /32)n 

(1  -  j82)2«-1 


YsMsNtp2* 

j= 0 


=  B(t,  rj,  <t,  r) 

► 

=  X(L,  6  »?) 


I  +  ( iesgnq )  ?? 

v'4I-(e  +  >;2) 


=  C(L,  f,  >?) 


£ J  _  eif[(n-2p+})A'  +  (n-m)i]e-im»( 


=  D(Y) 


(107) 


From  (61)  and  (107)  it  follows  that  the  final  partials  of  R  with  respect  to  the  Poincare 
variables  are 


dR  _  RdA  RdF  i?  &B  R  dX  RdC 

dL  ~  AdL+  fdL  +  V  B~dL  +  'x~dL  + 'q' C  dL 

dR  RdX  RdB  RdX  .  .RdC 

d^  =  J-dt+^BK+x~dt+^CK 

dR  _  RdF  RdX  RdC 

dg  ~  F  dg  +  'P]B  dg  +  X  dg  +  l9,C  dg 

dR  _  RdF_  RcFB 

da  f  da  ^  B  da 


\ 


>  (108) 


dR  _  RcW  RdB_ 

dr  ~  F  dr  +  ^  B  dr 

dR  _ 
dX  ~ 


RdD 
DdX  ) 
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From  the  preceding  pages,  we  find 


RdA 

A  dL 


=  -  (2n  +  2)  j 


ROF  _  R  (s')2  dF 
F  dL  f2c'  G  dc' 


RdF  R  (s')2  dF 
F  di  ~  F  ^  2  c'Gdc' 

RcW  -  *L,_  (3')2  dF 

F  drj  F  ^  2  c'  G  dc' 

RfrF  _  R  -a  dF 
F  da  F  4  c'  G  dc' 

R  dF  _  R  -t  dF 
F  dr  ~  f  4c'  G  dc' 


where 


dF 

dc' 


/  (n  -  r')!(n  +  r')! 

2”(n  —  m)\p\(n  —  p)\ 


* 


(x')  +  4(c'r 


d_K 

dx 1 


and 


P"^r>  (F)  obtained  with  (65)  and  (66), 


dP*  t'(-r') 


dx' 


obtained  with  (71)  and  (72). 


(109) 


(110) 


(HI) 
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RdB 

-M  R 

BdL~ 

2G 

RdB 

W\tR 

2  G 

RdB 

W\nR 

B  dn 

2  G 

RdB 

\v\R 

B  da  ~ 

2  \/GB 

.RdB 

\v\ (i  e  sgn  v 

B  dr  2 \/GB  ) 


Test  if  \u\  —  0.  If  yes,  put  all  (112)  to  zero  and  skip  (118b). 


R  OX  R  2 (G  -  L)  dX 
X  dL  -  x  (L  +  G)2  dp2 

R  dX  R  2£  dX 
J  di~  x  {L  +  G)2  dp2 

RdX  R  2?/  L  dX 

X  On  ~  x  (L  +  G)2  dp 2  , 


where 


dX 

dp2 


n  i  2  n  —  1 

(i  +  p2)  +  (T^) 


X  -  2(-l)kl 


(l+P2)"-1 

(1  -  p2)2”- 


7=0 


+  (-l)kl 


(i+iT 

(1  _/?2)2n-l 


^  MsNt  [j  p2  3  ~2\ 

7  =  1 


with  M.,,  IV(,  (Mv),  and  (IVhi)  given  by  (84)  and  (95),  or  where 


(112) 


(113) 

)  +  (Mv)N,\p23 

(114) 
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dX 
df3 2 


n 

,(i  +  /?2)  + 


2n  —  1 

(W5). 


X  +  2(-l)l?l 


(l+/^2)n~2 

(1  _/32)2n-l 


oo 


*  ^  [M(s  -  l,a  -  1;  b)M(i,7;  -6)  -  M(s,a;  6)M(i  -1,7-1;  -&))  /32j 

l=o 


+(-dw 


(l+l2)" 

(1-/S2)2"-1 


£M(i,Q;6)M(i,7;-6)[j/?2>-2] 

1=1 


(115) 


with  all  M  defined  by  (99)  and  (100). 


RdC  -\q\R 
lqlC  dL  ~  (L  +  G) 


M 


RdC 

c  di 


R  kl  L 

c  yj2(L  +  G)  [  yj'2 (L  +  G) 


> 


RdC_  _  R  \g | 

C  dt]  ~  C  y/2 {L  +  G ) 


(1  §)  + 


vC 

y/2  (L  +  G ) 


(116) 


If  |</|  =  0,  put  all  (116)  to  zero  and  skip  (118d). 


RdD 

Dd\' 


i  e(n  —  2p  +  q)R 


(117) 


All  ingredients  required  in  the  preceding  pages  are  defined  within  this  document;  they 
will  be  compiled  again  in  the  final  algorithm. 

Expressions  such  as  R  /  F  are  convenient  short-hand  notation.  Since  they  involve 
potential  zero  divisors,  they  should  be  computed  as  follows: 
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£  =  AD^XC^D 

F 

(118a) 

£  =  afbh~'xcmd 

skip  if  |^|  =  0 

(118 b) 

£  =  AFD^C^D 

X 

(118 c) 

£  =  AFB^XCM~lD 

skip  if  |</|  =  0 

( 118d) 

With  the  aid  of  (10S),  we  can  now  compute  the  six  partial  derivatives  which  we  will 

label 


OR 

*)« 


xj-L.i,  a ,  A',  ?/,  t 


(119) 


in  pc/ 


Compute 


0R\  V  (dR\ 

dii  )  v  Oij  ) 

J  '  "'ll  ,,=  ()  q—  —  oc  V  J  '  » 


y  =  L, 


(120) 


mpq 


D_ 

dy 


A  F-it 


+  complex  conjugate 


=  L. 


(121) 


Finally, 

Lit 

irr 

°2T 


dAFrr 

■  dAFiT 

A  ‘r  =  01 

d\' 

dAF2r 

dAFrr 

dy 

mi  ~  d£ 

dAFrr 

dAFrr 

dr 

Tir-  da 
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Combine  (122)  with  Lyddane’s  equation  of  motion.  The  later  may  be  found  in  Refer¬ 
ence  2,  page  A-26,  Fortran  lines  0648  to  0653. 


THE  ANALYTIC  SOLUTION 

PERTURBATIONS  IN  DELAUNAY  VARIABLES 
Each  Delaunay  variable  can  be  written  in  the  form 


L  —  Lq  4-  SL 


(123) 


where  SL  is  the  perturbation  due  to  the  tesseral  field.  The  Hamiltonian  is,  for  one  term 
in  n,  m,  p,  q , 


F  =  F0  -f  F%t 

=  +  (tH  +  7K)„mp,  (124) 


For  brevity,  the  subscript  n,  m,  p,  q  will  be  dropped  for  now.  The  familiar  canonical 
equations  of  motion  are 


L  =  F, 

1  =  -Fl 

> 

ii 

II 

1 

£ 

• 

(125) 

H'  =  Fh. 

cf 

1 

II 

•lie 

J 

We  take  the  time-derivative  of  (123),  apply  (125),  and  substitute  (124).  Since  L0  =  0. 
it  follows  that 
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6L  =  (7  Ri  +  7  Ri) 

6G  —  {jRg  +  7  Rg)  * 
6H'  =  (7/?*'  4-  7 Rw)  J 


(126) 


However,  for  /  we  find 

i  =  io  +  61 

=  -fl 

—  -Fol  -  {iRl  +7 Rl) 

=  -  [(Fol)[j0  +  Fqll6L\  -  (7 Rl  +7 Rl) 

Since 

io  =  -  (Fql)Lo 

there  follows 


6i  =  -{yRL+^Rr.)-Z^6L 

6g  =  -  (iRa  +  iRg) 

6h'  =  -  [~iRh'  +iRh') 


The  disturbing  function  may  be  written  in  the  form 


R  =  C(L.G,H')e 


1  v(l.gji') 


(12S) 


45 


NSWC  TR  89-109 


so  that 


Ri  =  ixp\R 


Substitution  into  (126)  yields 


which  integrates  to 


6L  =  iip\  (yR  —  jR) 


6L  =  Vi 
SG  =  Vs 


,~lR  +  •yR  \ 


ip' 

,  7  R  +  7  R 

ip1 


V' 


For  reasons  to  be  seen  later,  these  equatins  will  be  put  into  a  different  form, 
of  (128a),  we  write 


jR  +  jR-  — -  OyRi  -  7 Ri) 

V, 


(128a) 


(129) 


Because 


and  similar  equations  for  g  and  h' .  Then  (129)  becomes 
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In  order  to  solve  (127),  substitute  the  first  of  (129)  for  SL: 


This  integrates  to 


it  = 

iip'  L4  i{ip')2 


6g  =  - 


7  Rg  ~  iRg 


(131) 


fL,  iRh'-iRh' 

Oh  = - : - 


Above  equations  require  0j  and  0' .  Write  (28a)  in  the  form 


0'  =  e  j Ql  +  (n  —  2 p)g  +  mh'  —  em6t  +  (n  —  m)  — 


(132) 


There  follows 


0/  =  *Q 


(133) 


0'  =  e  \^Q(n0  +  h)  -f  (n  -  2p)g2  +  mh'2  -  emO 


(134) 


where 
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Upon  substitution  of  (133)  into  (131),  the  tesseral  perturbations  in  Delaunay  variables 
become 


6L  =  -4  (7 R,  -  7 Ri) 

xp* 

SG  =  j;(-rR!-lR,) 

6H'  =  ^-(7Rh.-*iHk,) 

61  =  J,  ^  +  i3Z7(Jp 

<5?  =  T7  (7-Rg  -  7-Rg) 

If)' 

6ti  =  -r -  (7^/7'  -  7#//') 

Ip* 


(135) 


PERTURBATIONS  IN  POINCARE  VARIABLES 


Expressing  the  tesseral  perturbations  in  terms  of  Poincare  variables  follows  a  sugges¬ 
tion  by  Cohen. 


Let 


X=(L,  G,  H\l ,  <7,  h'f 

and 


Y  =  ( L ,  a,  A',  r/,  r)T 
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Then  (135)  may  be  written  as 


X  =  X°  +  U 


( 1R  ~  7R\ 

V  ii>'  ) x 


(136) 


U  = 


h 

03 


Upon  making  the  canonical  transformation  to  Poincare  variables,  we  obtain 


Y  (X)  =  Y  (X°)  +  YxtU  f lR ,  .  7/?>) 

\  ixl>'  /  x 


=  Y°  +  YXtUYx 


,t/7R(Y)-7R(Y) 


iip'  (Y) 


iy0  ,  Tr,  JR-7R 


=  Y  +  U 


ii'1  Jy 


(137) 


where  the  temporary  tilde  emphasizes  that  R  and  ip1  are  now  functions  of  the  Poincare 
elements. 


Note  that 


0'(A ')  =  j'(Y) 


(13S) 
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and 


-  7/^ 


=  2  Im 


(139) 


namely  twice  the  imaginary  part  of  'yRy .  Substitution  of  (139)  into  (137)  yields  the  final 
expressions 


0 

A  L  —  —  Im(-yRy) 

ip1 


=  /m(7 R,,) 

ip* 


2 

A  a  =  —  Imi-yRr ) 

j/i' 


— ^  GcO  ti ^ 


Aii  =  Im(-yRt) 

0 


-2 

Ar  =  -A  /j/)(7/?.(T) 

ip1 


(HO) 


where  0'  is  obtained  by  evaluation  of  (134).  The  various  partial  derivatives  of  the  dis¬ 
turbing  function  7?  may  be  calculated  from  (108). 


PARTIAL  DERIVATIVES  OF  PSI  DOT 

The  prime  used  in  the  derivations  above  flagged  those  variables  which  were  “modified" 
Poincare  elements  or  functions  thereof.  This  prime  will  be  dropped  temporarily  so  that  it 
may  bo  available  for  a  different,  purpose. 
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As  seen  in  the  last  section,  we  may  write 


A  P  =  U 


2  Im(^R) 

4> 


(141) 


were  now 


PT  =  (L,  £,  cr,  A,  77,  t). 


Equation  (141)  yields 


A  P  =  U 


2  2 

—  Im(jRp) - —  Im('yR)il'p 

[xp  xp2 


(142) 


When  the  first  term  in  (142)  was  developed,  the  assumption  was  made  that  the  second 
term  is  negligible.  This  premise  proved  erroneous.  The  missing  terms  will  now  be  derived, 
following  Cohen’s  approach. 

Let  us  w'rite  (134)  in  the  form 

Xp  =  xp°  +  A  xp. 


where 


t/’°  =  eQn0  —  mO 


and 


A rj>  ~  -FlxTXpy. 


(143) 
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Here 


y  = 


I 

g 

eh 


„  h'CioRl 
F'  =  '  4 IW 


and 


0  =  e  Q/  +  (Q  -  ?)  g  +  em  (h  -  6)  +  (n  -  m)  ^ 


In  preparation  for  subsequent  steps,  we  make  a  canonical  transformation  to  the  vari 
ables 


=  (  L  ~  G  ]  ,  y  -  \  -g  -  eh\ 
\G-eHj  \  -eh  J 

and  obtain 


l*4Ci0R2E 

£  1 

3 

\(L-{L~G)-(G-tH)\2 

4  L6 

^5 
.  1 

i 

*-4 

[H  L-(L-G)  )  'J 

?//  =  e  Q\  +  q{—y  —  eh)  -f  (n  —  m  —  2 p)  (—eh)  —  emfi  +  (77,  —  777.)  — 


Ai>'  =  -eQF[  L  -  eqF[  L_G  -  evF[a_(l] 


(144 


where 


v  =  n  —  in  —  2  p. 
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The  leading  part  of  xj>p  is 


V’p  =  eQno^p, 


vhere  6  is  the  Kronecker  delta.  For  A  t/>p,  we  use 


Kypr  -  F[y  X'T  x'pt 


/I  0  0  0  0  0\ 

=  F[  ,  ,r  0  £  0  0  77  0 

11  \0  0  CT  0  0  T  J 


(145) 


For  the  sake  of  readability,  we  abbreviate 


L-G  =  a 


G-eH  =  f3 


Then 


I\L~3  (L  —  a)-3 


3 


(L-a-p 
V  L  -  a 


2 


1 


where 


A'  = 
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Then  the  required  partial  derivatives  become 

At/>(L,A)  =  -  (1, 0)  e  (QF'lLL  +  qF[L  a  +  vF'lL  ^  ' 

A =-(ff,r)c(QF,'M  +  qF'iaj  +  ^Fj'^)  ^ 


(146) 


After  several  pages  of  algebra,  the  six  partials  in  (146)  may  be  written 


6L-5  (L  -  cv)"5 


2(2 L-a)2  -L(L- a) 


3 


(^)2 


1 


-  6L~4  (L  -  a)~7  /3  [3  (51  -2  a){L  -  a  -  (3)  -  L{L  -  a)]  (147a) 


KFlL,a 


—3L~A  (L  —  a)~5  (7L  —  3a) 


3 


(^y 


+  6L~4  (L  -  a)"7  /3  [3  (4L  -  a)(L  -  a  -  f3)  -  L(L  -  a)\  (1476) 


=  18L-4  (L  -  a)-5  (21  -  a)(L  -  a  -  (3) 


+  6L~3  (L  -  a)~6  {L  ~  a  -2(3) 


(147c) 


12L-3  (L  —  a)~5 


-  6L~3  {L  -  a)~*  f3 


9 


(4^) 


(147d) 
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1F* 

jsrla,p 


—  18L~3  (L  —  a)-5 


(4H4) 


-6L~3(L-a)~<i(L-a-2ft) 


(147c) 


iF,'M  =  6L-3(L-a)-5 


(147f) 


where  A’,  a,  and  ft  were  defined  earlier.  However,  o  and  ft  may  also  be  computed  from 


«  =  5(f2  +  >/2)| 

/?  =  i(^  +  -2). 

With  the  abbreviations 

D\  =  e  (QF]L  ,  +  qFlL  a  +  vFu^) 
D2  =  e  (QF[L  o  +  qF'iQ  Q  +  vF[a > 
D\  =  (  ( QFxl,j  +  qFUxJ  +  vF\ti 

equations  (146)  become 

At/7,  =  ~£>i  Aj/»a  =  0 

Ai/’f  =  -££>2  At/’,,  =  ->/A>  1 

At/V  =  At/'r  =  —tDs 


(148) 


(149) 


(150) 


with  all  members  known. 
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Finally,  we  define 

2 

FAC  =  t* 

xl> 

FAD  =  FAC -Im^R) 

i' 

and  write  the  increment  to  Equations  (135)  in  vector  form: 

A  i'z 

—  A  rir, 

A  j/v 
-  A  rir 
A  rpL 
0 


RESONANCE 

CLASSIFICATION  BY  PERIOD 

Equation  (28a)  may  be  written  as 

=  (n  —  2p  +  q)l  +  (n  —  2 p)(g  +  h)  —  vh  -  emOt  -I-  (n  —  m)^,  (152) 

if  we  disregard  e  and  subscripts  for  the  time  being. 


Reintroducing  the  abbreviations 
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the  time  derivative  of  (123)  becomes 


ip  =  Q(n0  +  h)  -  m6  +  rg2  +  mh2 


(154) 


where  the  quantities  l2 ,  g2 ,  and  h2  are  the  first-order  secular  rates  of  these  elements,  due 
to  J2  =  —C2q  .  Formulae  for  the  rates  are  found  throughout  the  pertinent  literature  as, 
for  example,  in  Brouwer’s  (Reference  5)  celebrated  1959  paper: 


i2  =  \c{  1  -  3  cos2/)  y/\  -  e2  ' 

g2  =  ^C  ( 1  -  5 cos2  I) 
h2  —  C  cos  I 


(155) 


where 


—  1.5 J2R\ 

(1  -  e2)2  U 


Equation  (154)  is  used  to  examine  the  frequency  of  each  possible  combination  m,  Q.  r . 
While  the  period  of  ip  is  short  for  most  terms,  it  may  become  large  for  some  values  of 
Q  and  m.  When  the  period  is  long  ,  the  associated  perturbations  become  large,  and  the 
term  may  be  labeled  resonant.  To  be  more  specific,  a  term  is 


short-periodic  if 

9 

TS 

<H 

\ 

in  shallow  resonance  if 

9 

Td 

<H  < 

8 

TS 

> 

(15G) 

in  deep  resonance  if 

\ii'\  < 

9 

TO  ) 
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Theoretical  considerations  suggest  that  to  =  15  days  and  ts  =  0.5  days  are  reason¬ 
able  initial  trial  values.  However,  these  two  quantities  are  inputs  to  the  program,  to  be 
set  according  to  the  application  at  hand.  Since  the  inequalities  (156)  are  a  bit  hard  to 
visualize,  the  following  should  help  clarify  the  procedure:  if  the  critical  period  T ^ 


15"  <T 

the  term  is  integrated  numerically, 

0.5d  <T^  <  15d , 

the  term  is  evaluated  analytically, 

T+  <  0.5rf, 

the  term  is  ignored,  except  for  a 
few  cases  of  low  degree  and  order. 

Since  m  =  1  will  only  be  of  interest  in  dealing  with  geostationary  orbits,  we  will 
consider  all  m  from  2  to  nmax.  For  each  such  m  and  each  r  in  (—  nmax,  nmax) ,  obtain 
the  resonant  indices  Q  as  follows. 

Let 


a  = 


-  rg2 


n0  +  l2 


so  that 


Q  =  a  + 


no  +  l2 


Also,  let 


/?  = 


6 

n0  +  h 
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Consider  now  Figure  2. 


a 

_L 


a  +  — 
rs 


Q 


Figure  2.  Resonance  Conditions  for  Q 


If  there- is  an  integer  Q  in  the  interval  a  —  P/tq  to  a  +  P/tq ,  we  are  dealing 
with  deep  resonance.  If  there  are  any  other  integers  Q  in  (a  —  fi/rs,a  —  P/tq)  or  in 
(a  +  P/td,oi  +  P/ts ),  they  lead  to  shallow  resonance. 

Assuming,  once  again,  that  rs  =  0.5,  the  spread  in  possible  Qs  is  4 p.  If  4/?  <  1 , 
there  can  be,  at  most,  one  Q.  This  situation  may  occur  for  orbits  with  periods  of  less 
than  six  hours.  If  1  <  4/7  <  2,  there  can  be,  at  most,  two  Qs.  Such  conditions  may  arise 
for  orbits  with  periods  from  six  to  12  hours. 

We  have  now  identified  all  resonant  terms  with  indices  m,  Q,  r.  Since  the  algorithm 
operates  in  n,  m,  p,  q- space,  proceed  as  follows: 

n  assumes  all  values  from  nmin  to  nmax .  All 
m  have  been  identified;  they  must  also  obey  m  <  n. 
p  =  (n  —  r)  /2,  but  only  if  n  —  r  is  even,  and 
q  =  Q  -r. 

The  computer  will  write  two  files.  One  will  contain  all  n,  m,  p,  q  combinations  iden¬ 
tified  as  shallow  resonance,  the  other  all  those  that  are  in  deep  resonance. 


THE  AMPLITUDE  TEST 


Although  any  given  term  may  have  a  period  T ^  large  enough  to  classify  it  as  resonant, 
its  amplitutde  can  still  render  it  insignificant.  Since  a  first-order  theory  neglects  second 
order  perturbation  by  design,  it  serves  no  purpose  to  carry  terms  of  order  two  and  smaller. 
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Cohen  developed  a  test  in  which  the  estimated  amplitude  in  the  mean  anomaly  (or 
longitude)  is  compared  to  a  tolerance  of  order  two.  Considering  Equations  (135)  or  (140), 
he  approximates  the  amplitude  as 


Ampl 


(  J  AFX  \ 

V  Cn'4’  A 


7  A 


3  Qn0 


2n  ~b  2^ 


FX  +  slFsX  +  PlFX/,^ 


(157) 


where  the  disturbing  function  R  is  defined  in  (106)  and  (107),  and  where  CV  is  the 
normalization  coefficient 


C;v 


(tx  -  rn)!  (2 n  +  1)  (2  -  6om)  2 
(n  +  m)! 


A  is  defined  in  (61a),  F  is  the  normalized  inclination  function 


F  =  Cns^F 


with  the  latter  given  by  (44).  Moreover,  „Y  is  the  classical  Hansen  coefficient 


X  =  p^x, 


where 


0  = 


1  +  y/\  -  e2 


and  X  defined  by  (46).  The  two  products  of  partial  derivatives  in  (157)  are  found  to  be 


$lF,  - 


-2CnsM 


4  L-(C  +  if) 


.  ,  s2  dF 
\v\F--~ 
c  dc 


(158) 
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and 


PlXp 


4  L-i^+r]2) 


\q\X  +  2(3 2 


dX 

d(P 


(159) 


The  two  derivatives  dF/dc  and  dX/dfi2  have  been  defined  by  (76)  and  (96),  with 
(102)  as  an  alternative. 

If  the  amplitude  (157)  for  any  given  term  in  n,  m,  p,  q  exceeds  the  tolerance 


tol  = 


h  ( Re\ 2' 

2(1  -e2)2  U/  ’ 


(160) 


it  is  retained  as  a  resonant  perturbation  and  calculated  in  accord  with  its  classification.  If 
not,  the  term  is  dropped. 
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APPENDIX  A 
THE  FUNCTIONS  F  and  X 


A-l 
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The  literature  on  orbital  mechanics  contains  a  vast  collection  of  differenct  expressions 
for  the  inclination  function  and  the  Hansen  coefficient.  Some  of  the  classical  forms  of  these 
functions  are  presented  in  terms  of  summations  over  products  of  binomial  coefficients. 
Such  expressions  are  often  numerically  unstable,  and  many  significant  figures  are  lost 
when  calculations  are  made  for  large  values  of  n  and  m,  the  degree  and  order  of  a  tesseral 
term. 

In  order  to  overcome  this  difficulty,  and  to  reduce  execution  time,  a  variety  of  re¬ 
currence  relations  have  been  developed  and  published.  Cohen  and  Lyddane  perused  the 
literature  and  eventually  developed  serveral  methods  suitable  for  our  requirements.  Given 
below  are  the  presently  operational  versions  of  F  and  X .  Although  Lyddane  has  derived  ■ 
algorithms  more  elegant  than  our  current  X ,  we  will  defer  publication  until  numerical 
tests  are  complete. 


THE  INCLINATION  FUNCTION  F 


The  current  elegant  form  is  due  to  Cohen  who  credits  McClain  (Reference  G),  DiDo- 
nato,  and  Lyddane  with  the  inspiration.  Copying  Equation  (44), 


F'  ( c '  •' 

A  n  nip  v. 


,,  ,  (W  -  '■')!(»  +  *•)!  ,  ,,b  n0.6  ,  ,, 

)  =  (-!)  77777 - 77-77 - 77  c)  Pn_r,  (.r  ) 

2"(»  —  m)\p\(n  —  p)\ 


(Id) 


The  Jacobi  polynomials  P" ,br,  are  calculated  from  the  recurrence  relation  (45),  begin¬ 
ning  with  P'i  ’  as  function  of  P^l'b  and  P"'b .  Rather  than  derive  (1G1),  we  will  demonstrate 
that  this  expression  always  reduces  to  a  well-known  form  of  the  inclination  function.  Recall 
that 


r  =  n  —  2]>  b  —  -(n1  —  o) 


a  —  \w  —  ?’ 


h  =  |m  -f  r| 


,  1 


a  —  in  —  r  r  =  -(«  +  b) 


.r'  =  2(  r  )  -  1 


1  - 


A  -3 


+  t1 


4L-2(C  +  >r 
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The  Joacobi  polynomial  is  defined  by 


n  r>  +  b\s')2(n-r'~:i)(c')2] 


Now  let 


j  =  n  —  -(m  —  r  +  b)  —  k 


&2  =  k(j  =  0)  =  -(2n  +  v  ~  \^m  +  H) 


v  =  n  —  m  ~  2p, 


ki  =  k(j  =  n-  r')  =  +  \u\) 


Then 


pa.fc  _  Y'(_id(m-r-«)+^  n  +  \(a  b)  \/  n  \  (a  b)  \ 
n~r'  .  \k(m  -  r  +  a)  +  k)  \k(m  -  r  -  a)  +  k) 


*  (<j')2(fc-fci )(c')2n+l/-6-2* 


Substituting  this  in  FrJimp , 


pj  (  ,  1  (»-r')!(n  +  r,)i  ^  ■.  /  n  +  £(a  -  6) 

nmp( ,l  ’  2Il(n-?7i)!  p\(n  —  p)\  \J(m  -  r  +  a) + 


r  ( 77  i  —  ?■  +  a)  +  k 


*  (  n  ^(°  b)  )(s')2(k-k'\c')2n+,'-2k 
\i(77i  -  r  -  a)  +  kj 


Let  C;!mpk  be  the  coefficient  of  }2n+„-2k 
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We  will  substitute  for  a,  b  in  terms  of  m,  r  in  three  cases: 


(1)  r  >m.  Then 


a  =  r  —  m, 


b  =  r  +  m, 


l(a  +  6)=rr, 


-(a-b)  =  -m 


j  (-l)fc  (n  —  r)!(n  +  ?•)!  (n  -  m)!(n  +  to)! 

Cnmpfc  “  2  n(n  -  to)!  p!(n  -  p)!  A4(n  -  to  -  A-)!(m  ~r  +  k)!(n  +  r  -  fc)! 

(  —  l)*(n  +  to)!  /n  4-  r\  /  n  — r  \ 

2  np!(n  —  p)!  \  A:  )\n  —  m  —  k) 


(2)  m  >  r  >  —to.  Then 

a  —  m  —  r,  6  =  to  +  r,  -(a  +  6)  =  to, 

(  —  l)fc  (n  —  m)!(n  +  to)!  ^  n-r  \  f  n  +  r\ 

Lnmpk  =  2„(n  _  m),  p!(n-p)!  VTO-r  +  A:A  *  / 

(  — l)fc(n  +  m)!  /n  4-  r\  /  n-r  V 
2np!(n  —  p)!  \  k  J  \n  —  m  —  k ) 

(3)  —to  >  r.  Then 

a  =  to  —  r,  6  =  — ?n  -  r,  -(a  +  b)  —  -r,  -(a  -  b)  =  m 

J*  *-/ 


2  (°  ~  fe)  =  ~r 


j  (  —  1)*  (n  +  r)!(7r  —  r)!  4-  ?7i)!(n'—  m)\ 

nmpk  2 "(n  —  to)!  p!(?i  —  p)!  (?n  —  r  4-  k)\(n  4-  r  —  k)\k\(n  —  m  —  A-)! 

(  — l)fc(n  4- m)! /n  4- r\  /  n  — r  \ 

2"p!(n  -  p)!  \  k  J  \n  —  m  —  k) 
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In  all  three  cases,  C*  k  has  the  same  form  in  m,r.  Replacing  r  by  n  —  2p,  we  find 


(n  +  2n-2p\ 

2"p!(n-p)!^;  J  V  k  ) 


* 


(O 


2n+i/-2  k 


(162) 


where 


*i  =  M) 

k2  =  ^(2n  +  v  —  \2m  +  u\) 
s 

With  the  aid  of  Equation  (32),  (162)  is  seen  to  match  Equation  (12).  Q.E.D. 


THE  HANSEN  COEFFICIENT  X 


This  function  is  defined  by  Equation  (46).  Its  principal  components. 


u=0  x  7 


(163) 


u>=0  N  7 


(-by 


U’l 


(164) 


are  numerically  unattractive.  Lyddane  proposed  to  calculate  these  quantities  by  the  com¬ 
pact  recurrence  relations  shown  in  Equations  (49)  through  (51). 
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To  derive  Lyddane’s  expressions,  first  define  the  function 


Af(s,a;t)  =  (-1)*£ 

ti=0 


(165) 


Substitute  the  relation 


which  yields 


v=0  x  7 


=  M(s,a  —  1;  6)  —  M(s  —  l,a  —  1;6) 


which  is,  indeed,  Equation  (50).  Since  Nt  is  identical  in  form  to  Ms ,  except  for  the  signum 
of  b,  we  may  write 


Nt  =  M(t,r,-b). 


In  other  words,  the  relations  developed  for  M  can  be  used  to  evaluate  N  by  simply 
replacing  b  by  —b. 
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In  the  derivation  of  the  equations  of  motion,  derivatives  of  the  disturbing  function  are 
formed  which,  in  due  course,  require  differentiation  of  Ms  and  Nt  ■  With  the  aid  of  (165), 
we  write 


d_ 

db 


u=0  v  ' 


D=1  V 


bv-l 

iv  ~  !)• 


bv 


or 


^M(«,a;l)  =  -M(s-l,a-l;i)  (166) 

db 

In  a  similar  fashion  we  obtain 


~N(t,r,b)  =  M(t-  1,7  -  1;  -b)  (167) 

The  derivatives  (166)  and  (167)  are  subsequently  used  in  Equation  (102),  namely  the 
formula  for  dX  / d/?2 . 
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